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ABSTRACT 

A new method of infcrcncing of coupled stochastic nonlinear oscillators is described. The technique does not 
require extensive global optimization, provides optimal compensation for noise-induced errors and is robust in 
a broad range of dynamical models. We illustrate the main ideas of the technique by infcrcncing a model of 
five globally and locally coupled noisy oscillators. Specific modifications of the technique for infcrcncing hidden 
degrees of freedom of coupled nonlinear oscillators is discussed in the context of physiological applications. 
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1. INTRODUCTION 

Coupled oscillators are ubiquitous in nature. They are used to describe observed phenomena intensively over 
the years in many areas of science and technology including e.g. physics, 1,2 chemistry 3 and biology. 4 In this 
approach a complex system is characterized by projecting it onto a specific dynamical model of coupled nonlinear 
oscillators. However, there are no general methods to infer parameters of stochastic nonlinear models from the 
measured time-series data. Furthermore, in a great number of important problems the model is not usually 
known exactly from "first principles" and one is faced with a rather broad range of possible parametric models 
to consider. In addition, the experimental data can be extremely skewed, whereby important "hidden" features 
of a model such as coupling coefficients between the oscillators can be very difficult to extract due to the intricate 
interplay between noise and nonlincarity. 

As was pointed out by McSharry and co-authors, 5 deterministic inference techniques 6 consistently fail to 
yield accurate parameter estimates in the presence of noise. The problem becomes even more complicated when 
both measurement noise as well as intrinsic dynamical noise are present. 7 Various numerical schemes have 
been proposed recently to deal with different aspects of this inverse problem. 5, 7-12 A standard approach to 
this problem is often based on optimization of a certain cost function (a likelihood function) at the values of 
the model parameters that best reconstruct the measurements. It can be further generalized using a Bayesian 
formulation of the problem. 7, 9 Existing techniques usually employ numerical Monte Carlo techniques for 
complex optimization 11 or multidimensional integration 9 tasks. Inference results from noisy observations are 
shown to be very sensitive to the specific choice of the likelihood function. Consequently, the correct choice of 
the likelihood function is one of the central questions in the inference of continuous-time noise-driven dynamical 
models considered here. 

In this paper, we present an efficient technique of Bayesian inference of nonlinear noise-driven dynamical 
models from time-series data that avoids extensive numerical optimization. It also guarantees optimum compen- 
sation of noise-induced errors by invoking the likelihood function in the form of a path integral over the random 
trajectories of the dynamical system. The technique is verified on a broad range of dynamical models including 
system of five globally and locally coupled nonlinear oscillators. 

A specific example of inferencing stochastic nonlinear model from skewed time-series data is considered in the 
context of physiological research. In particular, we refer to the situation when the variability of the cardiovascular 
signals is modelled in terms of coupled nonlinear oscillators. ^' 13-15 At present there are no methods available 
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to infer parameters of the nonlinear coupling between oscillators directly from experimental time series data. 
Furthermore, in many situations it is important to perform such inference using univariate time series. This rises 
another important issue in nonlinear time-series analysis related to the inference of hidden dynamical variables. 
If a technique of inferencing of coupling parameters from hidden dynamical variables could be found it could 
provide new effective tool for estimation of the state of autonomous nervous control 16 and risk stratification of 
cardiovascular diseases. 17 The corresponding problem of inference of the coupling parameters of two nonlinear 
oscillators perturbed by noise from univariate time-series data will be considered in this paper. 

The paper is organized as follows. In the Sec. 2 the algorithm is introduced and its main features are 
compared with the results of earlier work. In the Sec. 3 the convergence of the algorithm is analyzed in the case 
of inference of coupled nonlinear stochastic oscillators. A modification of the algorithm that allows inference of 
hidden dynamical variables of two nonlinear coupled oscillators from univariate time-series data is considered in 
Sec. 4. 



2. THEORY OF NONLINEAR INFERENCE OF NOISE-DRIVEN DYNAMICAL 

SYSTEMS 

Consider iV-dimensional dynamical system described by set of nonlinear Langevin equations 

i(t) = f(x)+e(t) = f(x)+^(t), (1) 
where e(t) is an additive stationary white, Gaussian vector noise process 

<*(*)> = o, me(t'))=t)5(t-t% (2) 

characterized by diffusion matrix D. 

We assume that the trajectory x(t) of this system is observed at sequential time instants {tk\ k = 0, 1, . . . , K} 
and a series y = {yu = y{tk)} thus obtained is related to the (unknown) "true" system states X = {xk = x(tk)} 
through some conditional PDF p D (y\X). 

The most general approach to dynamical model inference is based on Bayesian framework (cf. 7 ). In the 
Bayesian model inference, two distinct PDFs are ascribed to the set of unknown model parameters: the prior 
Ppr(-M) and the posterior pp OS ^(M\y), respectively representing our state of knowledge about M before and 
after processing a block of data y. The prior acts as a regularize^ concentrating the parameter search to those 
regions of the model space favored by our expertise and any available auxiliary information. The two PDFs are 
related to each other via Bayes' theorem: 

Ppo S t(M\y)- Je{ylM)pMM)dM . (3) 

Here £(y\M), usually termed the likelihood, is the conditional PDF of the measurements y for a given choice 
M of the dynamical model. In practice, (3) can be applied iteratively using a sequence of data blocks y,y', 
etc. The posterior computed from block y serves as the prior for the next block y' , etc. For a sufficiently large 
number of observations, Pp 0S t(M\y , y' , . . .) is sharply peaked at a certain most probable model M = M* . 

The main efforts in the research on stochastic nonlinear dynamical inference are focused on constructing the 
likelihood function that compensates noise induced errors and on introducing efficient optimization algorithms 

( rf 5,7,9, 11). 

No closed form expression for the likelihood function that provides optimal compensation of the noise- induced 
errors was introduced so far for continuous systems. The ad hoc likelihood function 5 and their generalization to 
the conditional PDF for stochastic trajectories in maps 7,9 do not compensate the error in continuous systems, 
since they are missing the main compensating term (see below) . The problem of noise- induced errors in inference 
of continuous systems was considered in 11 and a general approach to constructing corresponding likelihood was 
outlined. However, the closed form expression for the likelihood that takes into account the leading compensating 
term was not found and instead an ad hoc expression for the likelihood function was used. 



A common draw back of earlier research is the use of extensive numerical optimization. This problem 
will become increasingly important when complex systems with the large number (hundred or more) of model 
coefficients are investigated. 

In the present paper we introduce a closed form of the likelihood function for continuous systems that provides 
optimal compensation for the noise-induced errors. We also suggest parametrization of the unknown vector field 
that reduces the problem of nonlinear dynamical inference to essentially linear one for a broad class of nonlinear 
dynamical systems. This allows us to write an efficient algorithm of Bayesian inference of nonlinear noise- 
driven dynamical models that avoids extensive numerical optimization and guarantees optimum compensation 
of noise-induced errors. 

In what follows in this section we describe the likelihood function, the parametrization, and the corresponding 
algorithm. 

2.1. The likelihood function 

It was pointed out in 11 the probability density functional for the nonlinear dynamical stochastic systems in 
general is not known. Instead one can use the probability density functional for random trajectories in such 
systems. We note that the path-integral approach has also proved to be useful in nonlinear filtration of random 
signals (see e.g. 18 )in the situations where standard spectral-correlation methods fail. 

Therefore we write the expression for the likelihood in the form of a path integral over the random trajectories 
of the system: 

i{y\M)= f x ' p (y\x)rM[x(t)]vx(t), (4) 

which relates the dynamical variables x(t) of the system (1) to the observations y(t). Here we choose U <C to < 
tic\t{ so that I does not depend on the particular initial and final states x(ii), x(if). The form of the probability 
functional Tm over the system trajectory x(i) is determined by the properties of the dynamical noise £(i). 19,20 

In the following we are focusing on the case of additive and stationary Gaussian white noise, as indicated 
in (1), (2). We consider a uniform sampling scheme tk = to + hk, h = {tx — to)/K and assume that for each 
trajectory component x n (t) the measurement error e is negligible compared with the fluctuations induced by the 
dynamical noise; that is, e 2 <C h(D 2 ) nn - Consequently, we use p (y\X) — > IlfcLo ^[y fe — x (*fc)] m (4)- Using 
results from 19 for ^vt[x(t)], the logarithm of the likelihood (4) takes the following form for sufficiently large K 
(small time step h): 

2 h K ~ 1 

--log£(y\M) = lndct D + - ]T [ tr *(y fe ; c) + (y fc - f (y fe ; c)) T D 1 y k - f (y fe ; c))] + N\n(2nh), (5) 

k=0 

here we introduce the "velocity" y^ and matrix 4>(x) 

Yfe =h~ 1 (y k+1 -y fe ), ($(x;c))„„/ = <9/„(x; c)/dx n >- 



It is the term tr«fr(yfc; c) that guarantees optimal compensation of the noise- induced errors in our technique 
and that distinguish our likelihood function from those introduced in earlier work. Formally this term appears 
in path integral as a Jacobian of transformation 21-23 from noise variables to dynamical ones. We emphasize, 
however, that this term is not a correction, but a leading term in inference as will be shown in the following 
sections. 

Note, that the optimization of the log-likelihood function (5) is in general essentially nonlinear problem that 
requires extensive numerical optimization. Below we introduce parametrization that allows to avoid this problem 
for a broad class of nonlinear dynamical models. In particular, a vast majority of the model examples considered 
in the earlier work on the nonlinear dynamical inference can be solved using this technique. Moreover, a large 
number of important practical applications can be treated using the same approach. 



2.2. Parametrization of the unknown vector field 

We parameterize this system in the following way. The nonlinear vector field f(x) is written in the form 

f(x)=U(x)c = f(x;c), 



(6) 



where U(x) is an N x M matrix of suitably chosen basis functions {£/„ m (x); n — 1 : N, m — 1 : M}, and c is 
an M-dimensional coefficient vector. 

The choice of the base functions is not restricted to polynomials, 0& (x) can be any suitable function. In 
general if we use B different base functions 0{,(x) to model the system (1) the matrix U will have the following 
block structure 



U = 



0i o 
0i 







02 
2 







0B 
0B 







(7) 



where we have B diagonal blocks of size N x N and M = B ■ N. 

An important feature of (6) for our subsequent development is that, while possibly highly nonlinear in x, 
f (x; c) is strictly linear in c. 

Eqs. (5) and (6) are two main ingredients that allow to solve problem of nonlinear stochastic dynamical 
inference analytically as shown in the following section. 

2.3. The algorithm 

The vector elements {c m } and the matrix elements {£>„„'} together constitute a set M = {c, D} of unknown 
parameters to be inferred from the measurements y. 

With the use of (6), substitution of the prior ppr(-M) and the likelihood £(y\M) into (3) yields the posterior 
Ppost(M\y) = const x exp{-S(M\y)], where 



S{M\y) = Sy(c,D) = ^y(D) - c T w y (D) + Vs y (D)c. 
Here, use was made of the definitions 

K-l 

p y (D) = h ifl V 1 Yk + K ln(det D), 

K-l 

Wy(D) = t^c w + hJ2 



U fc D y k — 



K-l 



H y (D) = Spr 1 + ft ^Tu^D-iUfc, 



(8) 

(9) 
(10) 
(11) 



fc=0 



where Ufc = U(yfc) and the components of vector v(x) are: 



N 



71=1 



8x n 



m = 1 : M. 



(12) 



The mean values of c and D in the posterior distribution give the best estimates for the model parameters 
for a given block of data y of length K and provide a global minimum to Sy (c, D). We handle this optimization 
problem in the following way. Assume for the moment that c is known in (8). Then the posterior distribution 



over D has a mean Dp 0gt = @y(c) that provides a minimum to Sy(c, D) with respect to D = D T . Its matrix 
elements are 

K-l 

(13) 



®f(c) ee 1 ]T [y* - U(y fc ) c] n [y k - U(y fe ) c 



fc=0 



Alternatively, assume next that D is known, and note from (8) that in this case the posterior distribution over 
c is Gaussian. Its covariance is given by 2y(D) and the mean Cp 0gt minimizes 5y(c, D) with respect to c 



c P ost = H y ( D ) w y( D )- 



(14) 



We repeat this two-step optimization procedure iteratively, starting from some prior values Cp r and Spr- 
it can be seen that the second term in the sum on the rhs of eq. (10) originating from tr$(y^) does not vanish 
at the dynamical system attractors (1), unlike the term (5) hJ2k=o D -1 y k corresponding to the generalized 
least square optimization. 24 Therefore both types of terms are required to optimally balance the effect of noise 
effect in {y k } (8) and provide the robust convergence. In the following section we analyze relative importance 
of both terms for the convergence of our algorithm. 

3. NUMERICAL EXAMPLES 

We verified the convergence and robustness of the algorithm on a broad range of dynamical systems. In this 
paper we will be specifically focused on the applications to the inference of coupled nonlinear oscillators. 

3.1. Five coupled oscillators 

Consider system of five locally and globally coupled van der Pol oscillators 



Xk = y k , 

5 5 

Vk =£fe(l -x k )yk -u k x k +y^r)kjXj + 7k,k-iXk-i(t) + jkM+iXk+iit) +yVfcj£j, (15) 



3 = 1 



We assume for simplicity that there is no observational noise and that the observed signal is y = (2/1, 2/2, 2/3, 2/4, J/s)- 
We note that for the model of coupled oscillators (15) parameters of the equations x k — 2/fe are known and do not 
have to be inferred. An example of a trajectory of (15) is shown in the figure 1(a) in projection on (x\, x 2 , £3) 
subspace of the configuration space of this system. We chose the following base functions 

<f>(l) = x 1 ; 0(2) =3/1; 0(3) =x 2 ; 0(4) = y 2 ; 0(5) = x 3 ; 0(6) = y 3 ; 
0(7) =x 4 ; 0(8) =2/4; 0(9) =x 5 ; 0(10) = y 5 ; 0(11) = x^; 0(12) = x 2 x 3 ; 
0(13) = X3X4; 0(14) = x 4 x 5 ; 0(15) = x 5 xi; 0(16) = xjyr, 0(17) = x\y 2 \ 
0(18) - x\y z - 0(19) - xj yi ; 0(20) = x 2 5 y 5 . 

Together with the elements of the diffusion matrix we have to infer 115 model coefficients. Example of the 
convergence of the coefficients of the 5 th oscillator to their correct values is shown in the Fig. 1(b). Results of 
the corresponding convergence for the 4 th oscillator are summarized in the Table 1. It can be seen from the 
Table that accuracy of estimation of the model parameters is better then 1%. Of a special interest for us is 
the compensation of the noise-induced errors. In the figure Fig. 2 we compare results of inference of one of the 
coefficients of the system (1) e\ for two different diffusion matrices D and 2D where matrix D was chosen at 
random 



D = 



0.0621 1.9171 0.4307 

0.5773 1.3597 0.3648 

1.9421 0.1099 0.1535 

1.9010 1.1997 0.0148 

0.4561 0.7863 1.5776 



0.0356 0.3113 

1.7559 0.3259 

0.7051 0.6268 

1.4443 0.0588 

1.9369 0.7153 



(16) 




Figure 1. (a) A projection of a trajectory of system (15) on (xi, xi, xz) subspace of its configuration space, (b) Conver- 
gence of the coefficients of the 5 th oscillator to the true values as a function of a number of blocks of data. We have 100 
blocks of data with 800 points in each block and the sampling time h = 0.02. oi = ei, 02 = — ei, £13 = — wi, 04 = ??i2, 
a& = f?i3, ae = f?i4, 07 = 7715, as = 715. The convergence of the five components of the diffusion matrix is shown in the 
insert. 



Table 1. Convergence of the coefficients of the 4 th oscillator of system (15). We have used 200 blocks of data with 
5000 points in each block. True values are shown in the first row, inferred values are shown in the second row, and 
corresponding standard deviations are shown in the last row. The accuracy of inference is within 5%. 



coefficients 


£4 


W4 




1142 


V43 


V44 


743 


745 


£>4l 


D42 


£43 


true value 


0.2 


-0.06 


-0.075 


0.24 


-0.23 


-0.2 


0.064 


0.095 


1.477 


2.316 


1.783 


inferred value 


0.199 


-0.062 


-0.069 


0.246 


-0.228 


-0.20 


0.066 


0.096 


1.477 


2.316 


1.782 


std 


0.005 


0.0033 


0.0046 


0.004 


0.0034 


0.0042 


0.001 


0.001 


0.001 


0.002 


0.002 



It can be shown that that without compensation term the estimator (14), (14) is reduced to the generalized 
least square (GLS) estimator. The Fig. 2 shows that the GLS estimator systematically overestimates the value 
of £1 and the larger is noise intensity the larger is the systematic error of the overestimation (see curves V 
and 2' for D and 2D correspondingly). By adding the term tr <f>(yk', c) we obtain optimal compensation of the 
noise-induced errors as shown by the curves 1 and 2 obtained for the same noise intensities. To see the effect of 
the compensation analytically it is instructive to rewrite the sum in the eq. (10) in the integral form 

,x(T) , ,T 

w y (D) = S pr 1 c p r+ / ij(y(t)fi>- 1 dy-- v(y h )dt, (17) 

Jx(t Q ) L Jto 

It can be seen from eq. (17) that for the attractor localized in the phase space the first integral is finite, since 
initial and final points of integration belong to the attractor. The second integral is growing when the total time 
of inference is growing. 

In particular, for a point attractor this integral is identically zero and the whole inference is due to the 
compensating term i f^v(y k )dt. This result is intuitively clear, since for the point attractor in the absence of 
noise the system will stay forever in the same point and no inference can be done. It is only noise that forces 
the system to move about in the phase space and makes it possible to perform inference. 
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Figure 2. Results of inference of the ei that were performed according to eqs. (10) - (14) (curves 1 and 2) are compared 
with the results of inference without compensating term tr €>(yfc; c) (curves 1' and 2') for different noise matrices d: d = D 
for (1) and (1'); d = 2 * D for (2) and (2') where D is given in eq. (16). 



4. INFERENCE OF TWO COUPLED OSCILLATORS FROM UNIVARIATE 

TIME-SERIES DATA 

As we have mentioned in the introduction in many real experimental situations the model is not usually known 
exactly from "first principles" and in addition, the experimental data can be extremely skewed, whereby impor- 
tant "hidden" features of a model such as coupling coefficients can be very difficult to extract due to the intricate 
interplay between noise and nonlinearity. 

A specific example of such experimental situation is inference of the strength, directionality and a degree 
of randomness of the cardiorespiratory interaction from the blood pressure signal. Such inference can provide 
valuable diagnostic information about the responses of the autonomous nervous system. 16 ' 25 However, it is 
inherently difficult to dissociate a specific response from the rest of the cardiovascular interactions and the 
mechanical properties of the cardiovascular system in the intact organism. 26 Therefore a number of numerical 
techniques were introduced to address this problem using e.g. linear approximations, 27 or semi-quantitative 
estimations of either the strength of some of the nonlinear terms 28 or the directionality of coupling. 29 ' 30 But 
the problem remains wide open because of the complexity and nonlinearity of the cardiovascular interactions. 

It is important to notice that simultaneous measurements of the cardiovascular signals is performed in different 
parts of the system (see e.g. 14 ). As a consequence the nonlinear characteristics of the oscillations are substantially 
modified in each signal and inference of nonlinear coupling parameters has to be performed preferably using 
univariate data e.g. blood pressure or blood flow signal only. The necessity to use univariate data in general 
poses serious limitations on the techniques of reconstruction and the problem become essentially nontrivial even 
in quasi-linear noise- free limit. 31-33 

In this section we investigate the possibility of extending our technique of reconstruction of coupled nonlinear 
stochastic oscillators to encompass the case of inference from the univariate time-series data in the context of 
physiological research. We note that this is a particular example of inference of hidden dynamical variables, 
which will be addressed elsewhere. An example of the actual signal of the central venous blood pressure (record 
24 of the MGH/MF Waveform Database available at www.physionet.org). The main features of the blood signal 
data is the presence of the two oscillatory components at frequencies approximately f r — 0.2 Hz and f c = 1.7 
Hz corresponding to the respiratory and cardiac oscillations. It can also be clearly seen from the spectra that 
the nonlinear terms including terms of nonlinear interaction (and cardiorespiratory interaction in particular) 
are very strong in this sample. We note that the relative intensity and position of the cardiac and respiratory 
components vary strongly from sample to sample with average frequency of the respiration being around 0.3 Hz 
and of the heart beat being around 1.1 Hz. To infer coupling parameters from the univariate blood pressure 




Figure 3. Example of the blood pressure signal (a) and of its spectrum (b) taken from the record 24 of the MGH/MF 
Waveform Database available at www.physionet.org. 



signal an important simplifying assumption can be used. Namely it is assumed that the blood pressure signal can 
be represented as the sum of the oscillatory components with the main contributions coming from the oscillations 
of the respiration and heart. 1 Accordingly we chose our surrogate data as a sum of coordinates of two coupled 
van der Pol oscillators s(t) = X\{t) + x 2 (t). It can be seen that the spectrum of s(t) (Fig. 4 (c)) reproduces 
mentioned above main features of the real blood pressure signal. 



xi = 2/1, 



Vi = - x\)yi 



2 2 
Ulfxi +OiiX 2 + ^ f3l,ijXiXj + ^ Tl.ij 
i,j=l '•?=} 



Xi]Jj 



(18) 



3=1 



Xl = 2/2, 



2/2 



£2(1 - 2:2)2/2 - uj%x 2 + a 2 x 2 + fo,ij x iXj + ^O x iVi (19) 



i,3 = l 



3 = 1 



Here noise matrix a mixes zero-mean white Gaussian noises and is related to the diffusion matrix D = a-a T . 

To infer parameters of nonlinear coupling between cardiac and respiratory oscillations we decompose "mea- 
sured" signal s(t) on two oscillatory components using a combination of low- and high-pass Butterworth filters 
representing observations of mechanical cardiac and respiratory degrees of freedom on a discrete time grid with 
step h = 0.02 sec. Obtained in this way signals So(t) and s\(t) are shown in the Fig. 4 (a) and (b) respectively.* 
To make this numerical experiment even more realistic the input signal s(t) was filtered before decomposition 
(using high-pass Butterworth filter of the 2 nd order with cut-off frequency 0.0025 Hz) to model standard proce- 
dure of de-trending, which is used in time-series analysis of the cardiovascular signals to remove low frequency 
non-stationary trends. We now use standard embedding procedure to introduce an auxiliary two-dimensional 
dynamical system whose trajectory z(i) = (zo(t), z\{t)) is related to the observations {s(tk)} as follows 



z n {tk) = 



S n (tk + h) - S n (t k ) 



where n = 1,2. The corresponding simplified model of the nonlinear interaction between the cardiac and 
respiratory limit cycles has the form (cf. with 15 ) 



'6,n*n-l 



h,r, 



+ bg^SnSn-1 + &10,nS„Z„ + bn,nS n Z n -l + &12,nS„_iZ„ + 6l3,nSn-l^n-l + bl4, n Z n Z n -l 
+ bl 5 , n slzl + he,nS n -lzl + £n{t), Tl = 0, 1. 



(20) 



where ^„ (t) is a Gaussian white noise with correlation matrix Q n n i . We emphasize that a number of important 
parameters of the decomposition of the original signal (including the bandwidth, the order of the filters) have to 

'Note the difference between actual oscillatory components of the original signal s(t): xi(t) and X2(i) and the compo- 
nents obtained using spectral decomposition with Butterworth filters: so(t) and si(t). 




O 1 2 f Hz O 1 2 f, Hz O 1 2 f, Hz 



Figure 4. Comparison of the power spectra of the inferred z(t) components of the signal (gray lines) with the original 
signal s(t) (black lines): (a) a low- frequency component of the signal So(t) obtained using low-pass Butterworth filter of 
the 5 th order with cut-off frequency 0.55 Hz (black line) is compared with the inferred signal zo(t)(grayline); (b) a high- 
frequency component of the signal si(t) obtained high-pass Butterworth filter of the 4 th order with cut-off frequency 0.55 
Hz (black line) is compared with the inferred signal z 1 (t) (gray line); (c) spectrum of the original signal s(t) = xi(t)+X2(t) 
(black line) is compared with the spectrum of the inferred signal z(t) = zo(t) + zi(t). 

be selected to minimize the cost (8) and provide the best fit to the measured time series {s(ifc)}. The parameters 
of the model (20) can now be inferred directly from the univariate "measured" time series data s(t). The 
comparison between the time series of the inferred and actual cardiac oscillations is shown in Fig. 4 and Fig. 
5. The comparison of the inferred parameters with their actual values is summarized in the Tables 2. It can 
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Figure 5. (a) Comparison of the inferred signal x(t) — x\(t)+x(2) (black solid line) with the original signal (black dotted 
line). In the insert fragments of both signals are compared with better resolution in time, (b) Comparison of the inferred 
phase space trajectory (x(t),y(t)) (black solid line) with the original one (black dotted line). To facilitate the comparison 
we have used the same initial conditions to generate phase space trajectory with exact parameters of the system and with 
inferred parameters of the system. 



Table 2. Comparison of the non-zero parameters of the second equation of the model (20) inferred from the univariate 
time series data s(t) with their actual values. We have used 1 block of data with 120000 points. True values are shown 
in the first row, inferred values are shown in the second row. 



coefficients 


62,1 


&2,2 


62,3 


b2A 


62,5 


&2,6 


h,9 


&2,14 


&2,16 




true value 


0.05 


-45.0 


-0.19 


0.25 





2.55 


0.2 


0.11 


-0.25 


0.2 


inferred value 


0.014 


-44.73 


-0.071 


0.17 


-0.081 


1.25 


0.415 


0.14 


-0.251 


0.17 



be seen from the Table that the inferred parameters give correct order of the magnitude for the actual values. 



The inferred values can be further corrected taken into account attenuation of the filters at different frequencies. 
We emphasize, however, that the technique of spectral decomposition of the "measured" signal is in principal 
non-unique. Moreover, in the actual experimental situation the dynamics of the physiological oscillations is 
unknown and can be only very approximately modelled by the system of coupled oscillators. Furthermore, the 
only criterion for the goodness of the spectral decomposition is the coincidence of the original and inferred signal 
and spectrum. For these reasons the estimation of the model parameters with the accuracy better then the 
order of magnitude does not improve the quality of the inferred information as will be discussed in more details 
elsewhere. 

In conclusion, we suggested new technique of inference of parameters nonlinear stochastic dynamical system. 
The technique does not require extensive global optimization, provides optimal compensation for noise-induced 
errors and is robust in a broad range of dynamical models. We illustrate the main ideas of the technique by 
inferencing 115 model coefficients of five globally and locally coupled noisy oscillators within accuracy 1%. It 
is demonstrated that our technique can be readily extended to solve selected problems of nonlinear stochastic 
inference of hidden dynamical variables in the context of the physiological modelling. We show in particular 
that the method allows one to estimate correct order of the magnitude of nonlinear coupling of two stochastic 
oscillators from univariate time series data. The framework of nonlinear Bayesian inference outlined in this paper 
can be further generalized to include errors of measurements and to solve problem of global inference hidden 
dynamical variables. 
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